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STRUCTURE FORMATION BY FIFTH FORCE: POWER SPECTRUM FROM N-BODY SIMULATIONS 
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ABSTRACT 

We lay out the framework to numerically study nonlinear structure formation in the context of scalar-field- 
coupled cold dark matter models (t^CDM models) where the scalar field (p serves as dynamical dark energy. 
Adopting parameters for the scalar field which leave negligible effects on the CMB spectrum, we generate the 
initial conditions for our Af -body simulations. The simulations follow the spatial distributions of dark matter and 
the scalar field, solving their equations of motion using a multilevel adaptive grid technique. We show that the 
spatial configuration of the scalar field depends sensitively on the local density field. The ipCDM model differs 
from standard ACDM at small scales with observable modifications of, e.g., the mass function of halos as well as 
the matter power spectrum. Nevertheless, the predictions of both models for the Hubble expansion and the CMB 
spectrum are virtually indistinguishable. Hence, galaxy cluster counts and weak lensing observations, which probe 
structure formation at small scales, are needed to falsify this class of models. 

Subject headings: cosmology: theory — methods: A^-body simulations 



1. INTRODUCTION 

The origin and nature of dark energy (Copeland et al. 2006) 
is one of the most difficult challenges facing physicists and cos- 
mologists at the present time. Among all the proposed models 
to tackle this problem, the introduction of a scalar field is per- 
haps the most popular. The scalar field, denoted by tp, should 
have no coupling to normal matter to be consistent with strin- 
gent constraints from experiments (Will 2006, and references 
therein), but could couple to the dark matter, therefore pro- 
ducing a fifth force between dark matter particles. This idea 
has gained a lot of interest in recent years because dark mat- 
ter physics are unknown, and such a coupling could alleviate 
the coincidence problem of dark energy (e.g., Amendola 2000; 
Chiba 2001; Chimento et al. 2003). Furthermore, it is com- 
monly predicted by low energy effective theories derived from 
a more fundamental theory. A specific and interesting possi- 
bility is the chameleon mechanism (Khoury & Weltman 2004; 
Mota & Shaw 2006), by virtue of which the scalar field ac- 
quires a large mass in high density regions and thus the fifth 
force becomes undetectable on short ranges, thus also evading 
constraints from the large-scale cosmic microwave background 
(CMB). Indeed, at the linear perturbation level, there have been 
a lot of studies about the coupled scalar field and f{R) gravity 
models (e. g., Li & Barrow 2007; Hu & Sawicki 2007). 

Nevertheless, little is known about these models on nonlinear 
scales. It is well known that the matter distribution at late times, 
i.e. z < 2 for cluster scales, evolves in a nonlinear way, mak- 
ing the behavior of the scalar field more complex and the linear 
analysis insufficient to produce accurate results that can be con- 
fronted with observations. For the latter purpose, the best way 
forward is to perform full A/-body simulations (Bertschinger 
1998) to evolve the individual particles step by step. 

A/-body simulations including scalar fields and related mod- 
els have been performed before (Linder & Jenkins 2003; Main- 
ini et al. 2003; Maccio et. al. 2004; Springel & Farrar 2007; 
Kesden & Kamionkowski 2006a, 2006b; Farrar & Rosen 2007; 
Baldi et al. 2008; Oyaizu 2008; Keselman et al. 2009; Li & 
Zhao 2009). For example, in the work of Maccio et al. (2004) 



the simulations included several effects due to the coupling be- 
tween dark energy and dark matter (e.g. modified gravitational 
constant, an extra dragging term in Newton's equations and 
time variable dark matter particle masses), but did not consider 
a spatial variation of the dark energy scalar field. The more 
complete simulation of the scalar field by Li & Zhao (2009) 
shows that this approximation is only good for a limited choice 
of parameters and the scalar field potential. Here we extend the 
work of Li & Zhao (2009). 

This paper is organized as follows: In § [2] we shall briefly 
review the general equations of motion for the coupled scalar 
field model introduced in Li & Zhao (2009), and present our 
specific choices of the coupling function and the scalar field 
potential. In §[3] we describe the formulae and the algorithm of 
the A/-body simulation, analyze the results of our coupled scalar 
field A/-body simulations, compare it with that of the standard 
ACDM model, and explain the physical origin of the new fea- 
tures. Finally, we conclude and discuss observational implica- 
tions in §|U 

2. THE COUPLED SCALAR FIELD MODEL 

2.1. The Model 

All properties of our coupled scalar field model can be de- 
rived from minimizing the action associated with the following 
Lagrangian density (the index a runs from to 3): 
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which includes the Ricci scalar R, and a dimensionless scalar 
field (p with a kinetic and an effective potential term. The latter 
is given by 

Veff(tp) = V(ip) - K(ip)C C DM (2) 

where the potential and the coupling function K(ip) are con- 
trolled by two dimensionless parameters, /i and 7, respectively. 
More rigorously, the potential V(ip) is 

(3) 
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and the coupling function n(tp) = 87rGexp(7</?), as given in Li 

6 Zhao (2009), where Ao is a constant on the order of the cos- 
mological constant, and G is Newton's constant of gravitation. 
Considering the non-relativistic, weak field limit of Eq. (0, 

VeffO) « A ^ + 8ttG(1 + 7^)p C DM, ~ (4) 
the meaning of this particular parameterization can be under- 
stood as follows: As the scalar field tp tends to minimize the ef- 
fective potential, the potential term Aop' 11 and the coupling (1 + 
jtp) to the CDM density (pcdm ~ ~£cdm in the non-relativistic, 
weak-field limit) lead to competing effects, favoring smaller 
and larger values of p, respectively. Q The balance of these two 
effects, minimizing the effective potential V e s, is controlled by 
the two dimensionless parameters /i and 7: fi is very small and 
controls the time when the effect of the scalar field (mainly ex- 
erting the finite-ranged fifth force on dark matter particles on 
galaxy cluster scales) becomes important for cosmology while 

7 determines how large it will ultimately be (Li & Zhao 2009). 
More specifically, the scalar field equation of motion is 



dV(tp) 

Up + — — + pcDM87rG7exp(7^) = 0. 

dtp 

Einstein's equations can be expressed as 

— -Gab = exp(7</?) / 0cDM«flMi + ?X 

ottU 
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where G„b is the Einstein tensor, and the right hand side is the 
energy-momentum tensor of the the scalar field and CDM with 
a four-velocity u a ; the scalar field's is given by 

"1 



UnGT* = V a pV b p- 
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Note that the energy-momentum tensors for the scalar field ip 
and the dark matter are not individually conserved due to their 
coupling, whereas their sum is. 

Eqs. ([5] © summarize all the physics that will be used in 
our analysis. An immediate application is the prediction of a 
uniform Hubble expansion. The model's expansion is com- 
pletely indistinguishable from ACDM for values of 7 ~ 0(1) 
and [i <C 1; the actual difference is on the order of O(p). Basi- 
cally, this is due to the large enough scalar's mass, forcing the 
field near the potential minimum, which itself is almost time- 
independent for /1 <C 1 ■ A quantitative explanation is given in 
Li & Zhao (2009). We now proceed to break the degeneracy 
via nonlinear clustering. 

2.2. The Nonrelativistic Equations 

The first step towards a numerical simulation is to simplify 
the relevant equations of motion in the non-relativistic and quasi- 
static limit (in the sense that the time derivatives can be safely 
neglected compared with the spatial derivatives). 

Li & Zhao (2009) showed that the scalar equation of motion, 
Eq. ((5]) and the Poisson equation can be simplified as 

pa 8^G 7 [pcdm-Pcdm]-Mo [p-^ 1 -^ 1 - 1 ] (8) 

ar 



» 4ttG [pcdm - Pcdm] - Ao [tp M - p p ] 



(9) 

Note that the above two equations have similar source terms, 
partly from matter and partly from the scalar field0 

' The dark matter Lagrangian £cdm specifies the geodesic flow for many point- 
like particles of four-velocity u a and density pcDM 

2 The notation 9 2 = — = 9 2 + 9 2 + d\ is defined with respect to the co- 
moving coordinate x such that V x = «V r , where a is the usual scale factor 
of the Universe. In the following, tp and pcdm denote the background val- 
ues of tp and pcDM: respectively. Although we have used the approximation 
Pcdm ewilP) ~ Pcdm f° r a simpler presentation, we keep the factor exp(7<p) 
as well as the potential given by Eq. {3} in the actual simulation. 




FIG. 1 . — Overdensity fields at z = for the (pCDM model with 7 = 1, fj, = 1(T 5 
(left) and the ACDM model (right). The former has developed more small-scale structure 
within the void. 



Finally, the equations of motion of the dark matter particles 
are also modified as 

dx p 

dt a 2 ' 



(10) 



dp 1 -± 
-T=— V x $-7V x </>, 
dt a 
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where the canonical momentum conjugate to the comoving co- 
ordinates x is p = a 2 x. Note that the two terms on the right hand 
side of Eq. (fTTT i correspond to gravity and fifth force, respec- 
tively (Li & Zhao 2009). The scalar field p is on the order of 
magnitude of /i, comparable to the dimensionless potential $. 
Eqs. ([8] |9j [10] [TTI ) are used in the code to evaluate the forces 
on the dark matter particles and to evolve their positions and 
momenta in time. 

The validity and limitation of the approximation present in 
the above equations, in particular neglecting the time deriva- 
tives, have been extensively discussed in Li & Zhao (2009). 
We emphasize that these approximations do not hold in linear 
regime where the scalar field's time dependence is essential for 
structure growth. However such terms have indeed been shown 
to be negligible on scales much smaller than the horizon scale 
(Li & Zhao 2009; Oyaizu 2008). To make our predictions more 
quantitative and rigorous compared to previous analyses (Mac- 
cio et. al. 2004; Kesden & Kamionkowski 2006a, 2006b; Farrar 

6 Rosen 2007), we now analyze the first A^-body simulations in 
the above framework. Considering the linear regime, Li & Zhao 
(2009) have already been able to constrain the parameters fi and 

7 to a fairly narrow range. Here we set 7 on the order of unity 
to force a significant ratio of the fifth force to gravity (~ 27), 
and explore the range 10~ 7 < [i < 10~ 5 , covering three orders 
of magnitude. Restricting ourselves to the above should suffice 
as the model is either essentially indistinguishable from ACDM 
or deviates too much from it (already at the linear level) beyond 
this parameter space, thus being of no further interest (Li & 
Zhao 2009). 

3. NONLINEAR STRUCTURE FORMATION 

In this section, we present some results of the first A-body 
runs and describe the qualitative behaviour of the coupled scalar 
field model. 

3.1. The N-Body Code 

We adapt the Multi-Level Adaptive Particle Mesh (MLAPM) 
code (Knebe et al. 2001) to include the scalar field, and its cou- 
pling to the dark matter A/-body particles. One benefit of the 
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adaptive scheme is that the majority of computing resources is 
dedicated to few high density regions to ensure higher resolu- 
tion, which is desirable since we expect the behaviour of the 
scalar field to be more complex there. 

The main modifications to the MLAPM code for our model 
are: 

1 . We have added a parallel solver for the scalar field based 
on Eq. ||8}. The solver uses a similar nonlinear Gauss- 
Seidel method (Briggs et al. 2000, Press et al. 1992) and 
the same criterion for convergence as the Poisson solver. 

2. The resulting value for (p of the first step is used to cal- 
culate the local mass density of the scalar field and thus 
the source term for Poisson's equation, which is solved 
using a fast Fourier transform to obtain the local gravita- 
tional potential $ [cf. Eq. ©]. 

3. The fifth force is obtained by differentiating ip, and the 
gravitational force is calculated by differentiating 4>, as 
inEqs. ([TOlIII). 

4. The momenta and positions of particles are then updated, 
taking into account both gravity and the fifth force, just 
as in normal A^-body codes. 

More technical details on the code, as well as how Eqs. ([HJ [9] 
[TOl ITTb are incorporated into MLAPM using its own internal 
units, have been given in Li & Zhao (2009) and will not be 
presented here. 

3.2. Numerical Results from the N-body Runs 

We have performed 6 runs of the modified code with param- 
eters 7 = 0.5, 1 and /i = 10~ 5 , 10~ 6 , 10~ 7 , respectively. For all 
these runs, there are 128 3 dark matter particles, and the simula- 
tion box size is chosen as B = 64/T 1 Mpc, with h being the usual 
dimensionless Hubble parameter and 128 domain grid cells in 
each direction. We assume a ACDM background cosmology 
which is a very good approximation for (Li & Zhao 

2009); in addition, we adopt present values for the fractional 
energy densities of dark matter and dark energy, JIcdm = 0.28 
and f^A = 0.72, and the normalization of the power spectrum is 
chosen as erg = 0.88. Note that the simulation does only take 
dark matter into account, baryons will be added in a forthcom- 
ing work to study the bias effect caused by the dark matter cou- 
pling. Given these parameters, the mass and spatial resolution 
of the simulation are 9.71 x 10 9 Mq and ~ 23.44/z" 1 kpc (for 
the most refined regions), respectively. This spatial resolution 
in high density regions is necessary and sufficent to precisely 
probe the scalar field in regions where the fifth force is consid- 
erably short-ranged. 

All simulations started at redshift z = 49. In principle, mod- 
ified initial conditions, i.e. the initial displacements and veloc- 
ities of particles which are obtained from a given linear matter 
power spectrum, need to be generated for the coupled scalar 
field model because the Zel'dovich approximation (Efstathiou 
et al. 1985) is also affected by the scalar field coupling. In 
practice, however, we find that the effect on the linear matter 
power spectrum at this high redshift is negligible (< ©(lO -4 )) 
for our choice of the parameters 7 and /i. Thus we simply use 
the ACDM initial displacements/velocities for the CDM par- 
ticles in our simulations, which are generated using GRAFIC 
(Bertschinger 1995), again using JIcdm = 0.28, tt\ = 0.72 and 
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FIG. 2. — Mass functions for 7 = 0.5 (upper panels) and 7=1 (lower panels) for 
different values of at z = and z = 1 . The ACDM mass function is also plotted as a 
(black) dot-dashed curve for comparison. 



o"8 = 0.88. An example of the final density field at redshift z = 
is shown in Fig.Q]for comparison with the ACDM simulation. 

We look for all virialized isolated haloes within our computa- 
tional volume using a Spherical Overdensity algorithm. For this 
purpose, we employ a time varying virial density contrast which 
is determined using the fitting formula presented in Mainini 
et al. (2003), and adopt the same virial density contrast for all 
models. In addition, we include all haloes with more than 200 
particles into the halo catalogue (see Maccio et al. 2008 for fur- 
ther details on our halo finding algorithm). Power spectra have 
been computed through a (fast) Fourier transform of the matter 
density field, computed on a regular grid Nq x Nc x Nc from the 
particle distribution via a Cloud-in-Cell algorithm (see Casarini 
et al. 2009). We set Nc = 256 which gives a maximum mode of 
k w 20/;Mpc _1 well above the simulation resolution. 

In Fig. |2] we show the mass functions for the runs with 
7 = 1.0,0.5 and fi = 10" 5 , 10" 6 , 10" 7 and the fiducial ACDM 
simulation at two output redshifts z = 1 and 0. The nonlinear 
matter power spectra of these models are displayed in Figs. [3] 
and[4] respectively. 

3.3. Interpretation of the Results 

The results of the A^-body simulations can be understood in- 
tuitively, as we shall discuss below. In general, a scalar field 
coupled to matter particles produces a fifth force [cf. Eq. ( fTTI il 
on the latter, which has a finite range m" 1 determined by the 
mass trip of the scalar field. If m v is small and almost constant 
across space then the fifth-force effect essentially leads to an 
increase in the effective gravitational constant which governs 
structure formation (Maccio et al. 2004). Li & Zhao (2009) 
have shown that for certain regions of parameter space and spe- 
cific choices of the potential, this is indeed a good approxima- 
tion. Mathematically this corresponds to neglecting the source 
terms starting with Ao in Eqs. d8l9t , hence the fifth force jVip 
is about a factor 27 2 times the gravitational force V$/a. 

In another situation, when the scalar field has a very steep 
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FlG. 3. — Ratios of calculated nonlinear matter power spectra for 7=1 and fi = 10 5 
(red), 10" 6 (blue) and 10" 7 (green) as well as for that of ACDM. Shown are results for two 
redshifts, z= 1,0. At large scales (small k) the curves converge to the horizontal curves 
(identical to 1, black dotted). Note that, using analytic results, the difference is expected 
to be small on both large and very small scales, and decreases at higher redshift. Error 
bars of future lensing observations are likely small enough to detect any deviation from 
ACDM on intermediate scales (k = 0.1- 10/jMpc" 1 ) at a 30% level. 



Du 



1 - 



Du 




-i 1 — 1 — 1 1 1 1 1 1 




k [h/Mpc] 

FIG. 4. — The same as Fig.[5] but for 7 = 0.5. 

potential, m v depends sensitively on the local matter density 
(Khoury & Weltman 2004) so that it almost resides at the min- 
imum of its effective potential 

V e ff(<p) = V(<A) + 87rGpcDMexp(7<^) (12) 

throughout space, i.e., <p ~ Ao/z/(87tGpcdm). This is known as 
the chameleon effect whose direct consequence is that in a high 
density environment, m 2 gets very heavy, 



: d 2 V eff /dif 2 



(13) 



and the fifth force becomes very short-ranged, with its effect be- 
ing suppressed due to 7V0 oc 7/uVpcbM- ^ n general the smaller 



/x is and/or the larger 7, pcdm are, the heavier becomes m v and 
thus the stronger the chameleon effect will be. Furthermore, 
since the value of ip inside a region also depends on its boundary 
condition, which in our case matches the background (p asymp- 
totically, we see that a smaller (p leads to a smaller <p and a 
heavier m v , and therefore to a stronger chameleon effect. 

There are several interesting features in Fig. f2] which can 
be understood schematically. First of all, our models produce 
more halos within the considered mass range than ACDM due 
to the enhancing effect of the fifth force. Secondly, a smaller 
[i means that the fifth force is more severely suppressed by 
the chameleon effect, and thus causes a small deviation from 
ACDM. Thirdly, a larger pcdm also means that the fifth force 
is more severely suppressed, and this is why at high redshifts 
the deviation from ACDM (for the same 7 and fi) is smaller. 
Fourthly, the influence of the parameter 7 is more complicated: 
A larger 7 will strengthen the chameleon effect, tending to sup- 
press the fifth force, but at the same time it increases the mag- 
nitude of the fifth force. In cases where chameleon effect is 
weak (e.g., /i = 10~ 5 ), however, we do see that a larger 7 leads 
to larger deviations from ACDM. 

Also note that the deviation from the ACDM mass function is 
more significant towards the low-mass end. To understand why 
this is the case, consider a mass range [Mo, Mo + AM]. At a cer- 
tain redshift, some halos which should have been in this range 
in ACDM indeed fall into the mass range (Mo + AM, 00) in our 
model as the fifth force accelerates the formation of structures 
(this tends to reduce the number of halos with mass > Mo as 
two halos which are separated in ACDM merge into one here), 
while some halos which should have been in the mass range 
< M in ACDM actually fall in the mass range [M ,M + AM] 
in our model (this increases the number of halos with mass 
> Mo). This effect is weaker for the largest halos because of 
competing effects due to merging of small halos. As the mass 
increases, the difference between the mass functions of the two 
models narrows down. 

In the matter power spectrum, we see something similar: 
Smaller /1 and larger pcdm (higher redshift) severely suppress 
the fifth force and lead to smaller deviations from ACDM; in- 
creasing 7 strengthens the fifth force, thereby causing large 
deviations from ACDM. Interestingly, the deviation becomes 
largest on intermediate scales: large scales are beyond the probe 
of the fifth force, and thus not significantly affected, while the 
density on small scales is high and the fifth force is suppressed. 

4. CONCLUSION 

We have presented a general framework to study nonlinear 
structure formation in coupled scalar field models, in particular 
the models of Li & Zhao (2009). While these models are virtu- 
ally indistinguishable from ACDM on both very large and very 
small scales, intermediate scales at low redshift (z < 1) relevant 
for galaxy clusters (~ 10 2 - 10 3 kpc) open a new window to test 
and constrain the interesting part of the parameter space. 

On these scales, the matter power spectrum is significantly 
increased compared to that of ACDM. Observationally, this would 
most likely appear as a change of erg on the order of 15-20% 
for models with 7 = 0.5- 1 and /i = 10~ 6 (see Fig. 2). Any 
variation of a% seems to be lower than 30% for current lensing 
measurements such as the CFHT Legacy Survey (e.g., Hoek- 
stra et al. 2006; Fig. 11 of Fu et al. 2008) over a rather lim- 
ited range; however, future surveys, such as the Kilo-Degree 
Survey (KIDS), will be able to measure the scale dependence 
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within the range k = 0.1 - 10/zMpc where the deviation of the 
models from ACDM is maximal. 
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